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We numerically study the symmetry breaking instabilities of vortex patterns in a rotating har- 
monic potential using a type of Ginzburg-Landau equation. The configurations of vortex lattices 
change markedly by the symmetry-breaking instabilities, and then, some vortices move away from 
the confinement potential, which leads to the annihilation of vortices. The symmetry-breaking 
instabilities and the instabilities of vortex nucleation determine the parameter region of stable vor- 
tex patterns. We verify that the symmetry-breaking instabilities also occur in a type of complex 
Ginzburg-Landau equation. 

PACS numbers: 03.75.Kk, 47.32.-y, 47.54.-r 

I. INTRODUCTION AND MODEL EQUATION 

Bose-Einstein condensates in trapped atomic eases at ultra-low temperatures have been intensively studied since 
their first experimental observation in 1995. G1US] ^ n ^ ne Bose-Einstein condensates, quantized vortices and vortex 
lattices were found experimentally in 1999. [Hll These vortex states have also been intensively studied theoretically. [f| 
0,| 

The dynamics of the Bose-Einstein condensates in a rotating potential is described fairly well by the Gross-Pitaevskii 
(GP) equation: 
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where m is the atomic mass, U = 1/kr 2 is the harmonic potential for the trapping, g = AirTi 2 a/m is a coupling constant 
characterized by the s-wave scattering length a, fl is the rotating frequency, and L z — —ih(xd y — yd x ). In this paper, 
we assume that h = 1 and m — 1 by performing a scale transformation, and we consider only a two-dimensional 
system. The vortex states can be stabilized by the rotation term — ilL z <f>. The solution with the lowest energy was 
approximately evaluated on the basis of the GP equation by Butts and Rokhsar. Q The linear perturbation S4> around 
a stationary solution </> obeys 
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If 8(f) is expressed as 8<p — {uje 
is obtained as 



)e and is substituted into eq. (2), the Bogoliubov-de Gennes equation 
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where j denotes the mode number and ujjis the eigenfrequency of the jth mode. The excitation of the surface mode 
is related to the nucleation of vortices. [1(| The excitation of the deformation mode of the vortex lattice is called the 
Tkachenko oscillation. [ll[ 

The energy and the total atomic number M = J |0| 2 <ir are conserved in the time evolution of the GP equation. 
Therefore, the lowest energy state is not attained by the direct numerical simulation of the GP equation from general 
initial conditions. Tsubota et al. performed a numerical simulation of the nucleation of vortices and the formation of 
a vortex lattice by introducing a damping term to the GP equation as (l2| 



-lv 2 +g\^\ 2 + U(v)- f i(t)-nL 2 



(4) 



where /i is the chemical potential and 7 is a damping parameter. The total atomic number M = J |^| 2 c?r is not 
conserved in this system, but they continuously adjusted the chemical potential /j,(t) so as to preserve the total 
number of condensates. 
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FIG. 1: Successive vortex patterns including (a) seven vortices, (b) six vortices, (c) five vortices, (d) four vortices, and (e) two 
vortices. 

In this paper, we study a type of Ginzburg-Landau equation to find stable vortex patterns. The model equation is 
written as 



dt 
dMx 
dt 



iv 2 - g\4>\ 2 - U{t) + m z + a(Mi(t) - M(t)) 
I3(M - M(t)), 



(5) 



where M(t) — J \<j>\ 2 dr, and Mo is a parameter representing the target total atomic number, Mi(t) is an additional 
variable, and a and (3 are additional parameters. The Ginzburg-Landau equation is obtained from eq. (4) by retaining 
only the damping term in the left-hand side of eq. (4) and setting the damping parameter to 7 = 1. The second 
equation for M\(t) acts as negative feedback to keep the total atomic number M(t) close to Mq- This model equation 
is an efficient model for obtaining stationary states with lower energy for a certain fixed value Mq of the atomic 
number, because the Ginzburg-Landau equation is a variational system. The Ginzburg-Landau-type equation eq. (5) 
can be directly obtained from the GP equation eq. (1) by a change of variable t — > it, in which t is interpreted 
as the imaginary time. However, in this paper, we consider eq. (5) to be a limit of eq. (4) for large 7 with the 
timescale changed from 7~ 1 i to t. In this interpretation, t is the actual time. An intermediate model between eq. (1) 
and eq. (5), which has a similar form to eq. (4), will be studied in §4. The stationary solutions to eq. (5) are the 
stationary solutions to the GP equation satisfying M(t) = Mq. These solutions are also stationary solutions to eq. (4) 
if the chemical potential \i is suitably adjusted. The chemical potential /j, is expressed as fi = a{M\(t) — Mq) for 
the stationary solutions to eq. (5). We have used this type of Ginzburg-Landau equation in a previous paper to 
find soliton solutions to the GP equation for attractive interaction. [l3T ] In numerical simulations in this paper, the 
parameter values of Mq = 210, a = 3, and j3 = 0.5 are used and a harmonic potential U(r) = l/2kr 2 with k = 0.25 
is assumed. We performed numerical simulations by the split-step Fourier method with 256 x 256 modes. 

The linear stability of the stationary solutions <p for eq. (5) can be studied using the linear equation for the 
perturbation 8(j): 
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if a{M\ — M) = /j, = const. Various types of linear instability can occur in the Ginzburg-Landau-type equation when 
the eigenvalue becomes positive. The positive eigenvalue corresponds to the situation where the eigenfrequency ujj 
is negative in eq. (3). When a surface mode exhibits an instability, new vortices are nucleated. When an instability 
occurs in a deformation mode of a vortex lattice, the configuration of the vortex lattice changes markedly. This is the 
main subject of this paper. 



II. SYMMETRY BREAKING IN POINT- VORTEX SYSTEMS 



Campbell and Ziff studied various vortex patterns in a rotating superfiuid confined in a cylinder of radius R using 
the energy of point vortices and their mirror images. [l4| They showed that many vortex patterns with various numbers 
of vortices are stable at the same frequency fl, and that multiple configurations are possible even for the same vortex 
number when it is larger than 9. That is, there are many vortex patterns, that are local minima of the energy We 
explain symmetry-breaking instabilities in the point vortex systems first, which will help with the understanding of 
the results of direct numerical simulations using eq. (5) described in §3. The energy / of a system of N point vortices 
in a rotating vessel with radius R and frequency Q, is expressed as, 
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where = y/ ' x 2 + y 2 is the distance between the center and the ith vortex located at (xi,yi), ri_j — 



yj (xi — Xj) 2 + (j/i — yj) 2 is the distance between the ith and jth vortices, r[ ■ — ^J(xi — R 2 Xj /r 2 ) 2 + (j/j — R 2 yj/r 2 ) 2 

is the distance between the ith vortex and the mirror image of the jth vortex located at (R 2 Xj/r 2 , R 2 yj/r 2 ), and 
is the external frequency. The equation of motion for the ith vortex is determined as 



dt dyi ' dt dxi ' 



To efficiently find stable stationary solutions satisfying dxi/dt = dyi/dt = 0, we study a variational system: 



dxi 
~~dt 
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This model equation corresponds to the Ginzburg-Landau-type equation eq. (5) for the GP equation. The direction 
of the velocity of a point vortex in the variational model is rotated by —tt/2 from that of the velocity given by eq. (8). 
We can study stable configurations of point vortices and their instabilities by numerical simulations using eq. (9). 

Firstly, we consider a single vortex located at (r, 0) . The position of the mirror image of the single vortex is 
(R 2 /r, 0). The energy / of the point vortex and the velocity v r are calculated as 



(The angular velocity vg obtained from eq. (8) is expressed as vg — l/(R 2 /r 2 — r) — Or.) The energy of the single 
vortex is / = InR at r = and / = — oo at r = R. The origin r = becomes a local minimum for O > l/R 2 . That is, 
a single vortex becomes a stable state for 51 > l/R 2 . This instability for < l/R 2 is a symmetry-breaking instability 
for a single vortex, because the most symmetrical position r = becomes an unstable position. We denote the critical 
values of il for the symmetry-breaking instability obtained from eq. (9) including N point vortices as fi c JV,o- ln our 
system confined in a harmonic potential U = (l/2)fcr 2 , the effective radius R can be estimated as R = (4Mo/7rfc) 1 / 4 
by the Thomas-Fermi approximation, where 4> becomes zero. In our model system with Mq = 210 and k = 0.25, the 
effective radius R and the critical value are evaluated as R = 5.71 and fid.o ~ 0.0306. In the time evolution of eq. (9), 
a single vortex stays at the center for fi > fi c i,o, but it becomes unstable and the vortex moves away from the center 
for O < Q c ifi. This instability is due to the force from the mirror image. 

Next we consider two vortices located at (ri,0) and (— r 2 ,0). The vortex position in the case of a symmetrical 
solution satisfying n = r 2 = d/2 is obtained as a solution of = 8d 2 /(16R 4 — d 4 ) + 2/d 2 , which is approximated to 
d = \J2/Q., if the mirror-image force is neglected. The stationary solution r\ = r 2 = d/2 disappears for f2 < 0.065 
when R = 5.71. However, a symmetry-breaking instability occurs at a larger critical value of f2 C 2,o = 0.085 in the 
time evolution of eq. (9). For £1 < 0.085, one vortex moves to the center and the other one moves towards the wall 
at r — R owing to the symmetry-breaking instability. The symmetry breaking implies that the energy / is not a 
local minimum at r± = r 2 = d/2 for < Q C 2,o- The symmetry-breaking instability does not occur if the image force 
is neglected. In the time evolution of eq. (9), the creation and annihilation of point vortices do not occur, and the 
distance r from the center of one vortex can become larger than R; in this case, eq. (9) has no physical meaning. 

Similar symmetry-breaking instabilities also take place for systems including a large number of vortices. Figure 1 
displays vortex patterns that are expected to appear after successive symmetry-breaking instabilities from the first 
pattern, which has 7 vortices. This sequence of vortex patterns was obtained by a numerical simulation using eq. (9) 
by gradually decreasing the parameter fi. In the actual numerical simulation, the variables Xj and yj for the jth 
vortex were removed when rj reached R. This is because the vortex is expected to disappear at the wall, where 
r = R. We have numerically confirmed that these successive symmetry-breaking instabilities do not occur if the 
forces from the mirror images are neglected. The critical values of the successive symmetry-breaking instabilities 
were numerically found to be Q C 7.o — 0.231, f2 C 6.o = 0.209, f2 C 5.o = 0.181, and f2 C 4,o = 0.143 for R = 5.71. We can 
evaluate another value of the effective radius, R = 7.5, which was found by a direct numerical simulation of eq. (5) 
as a radius satisfying \<f>(r)\ <~ for r > R. If the cylinder radius is assumed to be R — 7.5, the critical values of the 
successive symmetry-breaking instabilities are ri C 7 : o = 0.133, £l C 6,o = 0.116,f2 C 5,o = 0.1, fl C 4,o — 0.083, Q C 2,o = 0.049, 
and fl c i,o = 0.017. The critical values depend strongly on the radius R, because the mirror force is essential for the 
symmetry-breaking instabilities. The mirror symmetry is broken at the instability of vortex patterns (a), (b), (c), and 
(e). On the other hand, for the vortex pattern (d) with iV = 4, the square configuration becomes a rhombic pattern 
at the symmetry-breaking instability. After the instability, two vortices move outwards and disappear at r = R and 
two vortices survive, as shown in Fig. 1(e). 



/ = (1/2) \n(R 2 - r 2 ) + (l/2)ftr 2 , v r = l/(R 2 /r - r) - fir. 



(10) 
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FIG. 2: (a) Profile \<j>\ at the section y = L/2 for Q. = 0.14. (b) Time evolution of \(j>\ at the section y = L/2 for Q = 0.12. (c) 
Time evolution of the vortex position X v for Q, = 0.12. 




FIG. 3: (a) Distance d between two vortices for N = 2. (b) Distances d between two vortices for N = 4 and N — 7. 



III. SYMMETRY BREAKING AND VORTEX ANNIHILATION IN THE GINZBURG-LANDAU-TYPE 

EQUATION 

In this section, we study vortex patterns and their symmetry-breaking instabilities by numerical simulation using 
eq. (5) with g = 1. The system size is L x L = 30 x 30. The center of the harmonic potential is assumed to be 
(L/2, L/2). In contrast to the point-vortex system, our vortices have a finite core size. Our vortices are not confined 
by a cylindrical wall but by a harmonic potential. Also, the nucleation and annihilation of vortices occur naturally in 
this system in contrast to the point-vortex system in §2. 

We first performed the numerical simulation of a single vortex, which was initially set near the center. Figure 2(a) 
displays \4>\ for £1 — 0.14 in the cross section at y = L/2. A vortex is stably located at the center for the frequency 
£1 = 0.14. The modulus \4>\ is at the position of the vortex, and 4> is also almost for r > R ~ 7.5 owing to the 
confinement by the harmonic potential. Figure 2(b) displays the time evolution of \<f)(x,y)\ at the section y = L/2 
by eq. (5), and Fig. 2(c) displays the time evolution of the vortex position when £1 = 0.12. The single- vortex state 
is unstable at this frequency, and the vortex moves away from the center and finally disappears. We denote the 
critical values of $1 for the symmetry-breaking instability in eq. (5) including TV vortices as Q c n- The critical value 
is £l c i ~ 0.13, which is interpreted to correspond to fi c i,o = 0.0306 for R — 5.71 or Qci,o = 0.017 for R = 7.5 in 
the point-vortex theory. On the other hand, new vortices are nucleated from the instability of the surface mode at 
il = 0.31. We denote the critical values of for the vortex nucleation in eq. (5) including N vortices as £l n N- That 
is, fl n i — 0.31. 

Next, we consider a system including two vortices, which are initially located at (xi, L/2) [x\ > L/2) and (xz, L/2) 
(X2 < L/2). Figure 3(a) displays the distance d = x\ — X2 between the two vortices in the stationary solution as 
a function of Q obtained by a long time numerical simulation using eq. (5). The dashed curve is the theoretical 
curve estimated from f2 = 8d 2 /(16i? 4 — d A ) + 2/d 2 by the point-vortex theory. Here, R = 7.5 is used; however, the 
effect of the mirror image is very small and the curve is well approximated at d = y/2/Cl for the parameter range 
0.2 < ft < 0.31. In the numerical simulation of eq. (5), the symmetry-breaking instability occurs at f2 C 2 = 0.19, which 
corresponds to f2 C 2,o = 0.085 for R = 5.71 or fi C 2,o = 0.049 for R = 7.5 in the point-vortex theory. Although similar 
types of symmetry breaking occurs according to eq. (5) and eq. (9), both the distance between the two vortices and 
the critical value of the symmetry breaking are quantitatively different from those in the point-vortex theory. Figure 
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FIG. 4: (a) Time evolution of \tf>\ at the section y = L/2 for SI = 0.18. (b) Time evolution of the two vortex centers xi and X2 
for fi = 0.18. 




FIG. 5: (a) Vortex pattern at t — 100 including four vortices when SI = 0.21. (b) Vortex pattern at t — 400 when fl = 0.21. 
(c) Time evolutions of the x- and y-coordinates of the vortex positions. 



4(a) displays the time evolution of \(f>\ at y = L/2 for S7 = 0.18. Figure 4(b) displays the time evolutions of the 
x-coordinates of the left and right vortices. It is clearly seen that one of the vortices moves to the center and the 
other one moves outwards and disappears. A single-vortex state appears after a long time, which is located at the 
center. On the other hand, the vortex nucleation occurs at fl n 2 = 0.32 for this system with N = 2. 

We have studied a system including four vortices forming a square and a system including seven vortices forming 
a triangular lattice, whose configurations are respectively similar to Figs. 1(d) and 1(a). Figure 3(b) displays the 
distances between the nearest-neighbor vortices for the systems with N = 4 and N = 7. The dashed curves are 
theoretically estimated by the point-vortex theory, where R = 7.5 is used. The theoretically estimated values are 
smaller than the numerical results. The symmetry-breaking instability occurs at S1 C 4 = 0.23 for the system with 
N = 4 and at S7 C 7 = 0.29 for the system with N — 7. Vortex nucleation occurs at Sl„4 = 0.35 for the system with 
TV = 4 and at fi n 7 = 0.40 for the system with N — 7. The critical value £l n N of the vortex nucleation increases with 
the vortex number N. The vortex pattern with vortex number N is stable in the parameter region between fl C N and 
QnN • The parameter regions mutually overlap for different values of N. We can therefore conclude that multiple 
vortex patterns are also stable in our system for the same frequency SI from the relation bwteen the critical values 
fl C N and fl n N for various N. 
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FIG. 6: Vortex patterns at (a) t = 100, (b) t = 150, (c) t = 1000, and (d) t = 1750 for a system with N = 7 when fi = 0.28. 

(a) (b) 
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FIG. 7: (a) Trajectories of four vortices obtained from eq. (11) when fi = 0.081, R = 7.5 and j = 0.1. (b) Trajectories of four 
vortices by the complex Ginzburg-Landau-type model eq. (12) when Q = 0.21 and 7 = 0.1. Solid curves correspond to the 
vortices taht survive and dashed curves correspond to the vortices that move outwards and disappear. 



Figures 5(a) and 5(b) display snapshot vortex patterns at t = 100 and t = 400 as a result of the symmetry-breaking 
instability when = 0.21 for the system including four vortices initially located on the lines x = L/2 and y — L/2. 
The two vortices on the line y = L/2 move outwards, and the two vortices on the line x — L/2 survive, which is 
qualitatively consistent with the point-vortex theory. Figure 5(c) displays the time evolutions of the x-coordinates 
Xi and X 3 of the vortex positions on the line y = L/2 and the y-coordinates Y 2 and Y4 on the line x = L/2. The two 
vortices on the line x = L/2 approach each other and the other two vortices on the lines y — L/2 move away from 
each other. 

Figures 6(a)-6(d) display four snapshots of vortex patterns at (a) t — 100, (b) t — 150, (c) t = 1000, and (d) 
t = 1750 for a system with N = 7 when Q = 0.28. As a result of the symmetry breaking, one vortex moves away at 
t = 150, then the other two vortices also move away. Four vortices survive at t = 1000; the configuration of the four 
vortices changes markedly and they finally form a square pattern at t = 1750. 

These numerical results show that symmetry-breaking instabilities also occur according to eq. (5), which are qual- 
itatively the same as those found in §2 for the point vortices. However, the critical values £l C N,o m the point-vortex 
system given by eq. (9) and £I c n given by eq. (5) are rather different. This is partly because the core size of the 
vortex is finite in the system given by eq. (5), in contrast to the point vortices with infinitesimal core size. As the 
nonlinear parameter g is increased, the core size decreases as for large g from an argument based on scaling. 

That is, if the core size is denoted as £, (1/2)V 2 ~ l/(2£ 2 )</> ~ g\<t>\ 2 4> f° r large g in eq. (1) with % = m = 1. 
The point-vortex theory is therefore expected to give a better approximation in the case of large g. Furthermore, 
the harmonic potential is used in the system based on eq. (5) to confine vortex patterns in contrast to the circular 
confinement by a cylindrical wall with radius R for the system given by eq. (9). 
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IV. SYMMETRY BREAKING IN A COMPLEX GINZBURG-LANDAU-TYPE EQUATION 

The symmetry-breaking instability does not occur in the GP equation eq. (1) or in eq. (8) for the point vortices, 
because there is no energy dissipation. If some dissipation terms are included, the symmetry-breaking instability is 
expected to occur. We have verified this using intermediate models between eq. (8) and eq. (9), and between eq. (1) 
and eq. (5). A model equation for the point vortices with a small amount of dissipation is expressed as 

— = 7o— -71— — = -7o— - 7i— (11) 
dt dyi dxi ' dt dxi dyi ' 



where 70 = 1/yl + 7 2 , and 71 = 7/yl + 7 2 using a parameter 7 representing the degree of dissipation. This model 
equation has the form of a linear combination of eqs. (8) and (9). Similarly, we can construct a model equation from 
eqs. (1) and (5) as 



-qT = H70+71) 



iv 2 - g\cp\ 2 - U(v) + flL 2 



cj) + -f 1 a{M 1 {t)-M{t))(j}, 



^ = 0(Mo-M(t)), (12) 

where 70 = 1/ \J\ + 7 2 , and 71 = 7/ ^/l + 7 2 . This is a model equation with a similar form to eq. (4) and is interpreted 
as a type of complex Ginzburg-Landau equation. When 70 = 1 and 71 = 0, eq. (11) is reduced to eq. (8), and eq. (12) 
reduced to eq. (1) when h = m = 1. When 70 = and 71 = 1, eq. (11) is reduced to eq. (9) and eq. (12) reduced to 
eq. (5). We have performed a numerical simulation of vortex patterns with four vortices initially located on the lines 
of x = L/2 or y = L/2. Figure 7(a) displays the trajectories of the four vortices using the point-vortex model eq. (11) 
when fl = 0.081 obtained from R = 7.5, and 7 = 0.1. Figure 7(b) displays similar trajectories of the four vortices 
in the complex Ginzburg-Landau-type equation eq. (12) when fl = 0.21,7 = 0.1, and g = 1. (For comparison with 
Fig. 7(a), the origin of the coordinates is shifted to (0,0) in Fig. 7(b).) The symmetry-breaking instability occurs 
for the vortex patterns and two vortices move outward and disappear. The other two vortices move inwards, and 
find new stationary stable positions. These symmetry-breaking instabilities are essentially the same as those obtained 
from eqs. (5) and (9). However, the rotational motion (spiral-like motion) of vortices is observed in these models, 
although the motion is rather complicated. The rotational trajectories are characteristic of the motion of vortices 
given by the GP equation eq. (1) and eq. (8) for the point vortices. These numerical simulations suggest that the 
symmetry-breaking instabilities are not only found in the imaginary models eqs. (5) and (9), and the instabilities are 
realistic if some energy dissipation is involved in the system. Rotational or spiral-like trajectories seem to be more 
realistic than the straight trajectories shown in Figs. 4 and 5, and such trajectories might be found in experiments. 
The straight trajectories numerically obtained by eq. (5) and shown in Figs. 4 and 5 are expected to appear in the 
limit of strong dissipation. 

V. VORTEX PATTERNS NEAR THE LOWEST-LANDAU-LEVEL REGIME 

The GP equation eq. (1) with m = fi = 1 is expressed as 



M_1 



d .. V f.B 



i- — h Qx) + i- fly 

ay J \ ox 



(13) 



in the special case when g = and Q — \/k. There is a family of vortex solutions to eq. (13): 

4> = {x - x + i(y - y )} cxp[-fl/2{(x - x ) 2 + (y - yo) 2 }} exp[ifl(x y - yox)}, (14) 

where (xo,yo) is the core position of a single vortex. This is equivalent to the wave function for an electron in a 
magnetic field at the lowest Landau level. We therefore call the parameter region of g ~ and Q ~ y/k the lowest- 
Landau-Level regime (L-L-L regime). A large number of vortices are nucleated near £1 — \fk if g is not small. In the 
parameter region, Ho proposed a solution for a vortex pattern including TV vortices: 

4> ~ Ilf =1 [(z - Xj ) + i(y - %)] exp{-r!/2(x 2 + y 2 )}, (15) 

where (xj,yj) is the core position of the jth vortex. [l5| However, such a Gaussian behavior was not observed 
for the envelope of <j) with a large number of vortices. It is considered that the Thomas-Fermi form is a better 
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FIG. 8: (a) Profile \<$>(x,y)\ of a single-vortex pattern at y = L/2 for g = 0.005, Q = 0.49, and M = 210 using eq. (5). The 
dashed curve denotes 4.09|x — L/2\ exp{— 0.25(x — L/2) 2 }, although it almost overlaps with the solid curve of \<j}\. (b) Profile 
of \((>(x,L/2)\ of a two-vortex pattern for g = 0.005, = 0.49, and M = 210. The dashed curve denotes 1.63|(a; - L/2) 2 - 
1.47 2 |exp{-0.205(a;-L/2) 2 }, and the dotted curve denotes 1.63|(x - L/2) 2 - 1.47 2 | exp{-0.25(x - L/2) 2 }. (c) Comparison of 



contour plots of the phase of < 
are drawn. 



i obtained from eq. (5) and the ansatz eq. (16). The contour lines of (j> — n7r/4 with n = 0, 1, 2, • • • , 8 



approximation. [161 ] The Thomas-Fermi approximation is suitable when g is not small. We have studied vortex 
patterns with a few vortices in the L-L-L regime, because the wave function of a single vortex is explicitly given as 
eq. (14), and a theoretical approach different from the point-vortex theory might be possible. 

We have performed numerical simulations using eq. (5) by setting a few vortices as an initial condition for sufficiently 
small values of g at ft — 0.49, k = 0.25 a = 3, /3 = 5, and M = 210. Figure 8(a) displays the profile of \<fr(x, y)\ at the 
section y = L/2 of a single vortex when g = 0.005 and £1 = 0.49. The dashed curve denotes A\x — L/2| exp{— 0.25(a; — 
L/2) 2 }, where A is calculated as i/M/(4tt) - 4.09 fr om the normalization condition. The vortex solution given by 
eq. (14) is a good approximation. 

Figure 8(b) displays the profile of \cf>\ at the section of y — L/2 for a system including two vortices at g — 0.005. 
The distance between the two vortices is 2.93 at g = 0.005. The dotted curve is the approximation \<f>\ = A\(x — 
L/2) 2 — Xq\ exp{— 0.25(ie — L/2) 2 }, which is obtained from the ansatz by Ho: 

<j> = A{x - L/2 - .T + i(y - L/2)}{x - L/2 + x + i{y - L/2)} cxp[-0.25{(a; - L/2) 2 + (y - L/2) 2 }}, (16) 

and the vortex position xo — 1-47 obtained from the numerical simulation is used. The ansatz given by eq. (16) is 
not satisfactory even for the small value of g. The dashed curve, which almost overlaps with the numerical result, is 
\4>\ — A\(x — L/2) 2 — Xq\ exp{— 0.205(x — L/2) 2 }. This Gaussian approximation is rather good. Figure 8(c) displays 
a contour plot of the phase of <j) given by tan -1 (Re</>/Im</>). The solid curve denotes the numerical results and the 
dashed curve is obtained from the ansatz eq. (16). These figures imply that the ansatz is unsatisfactory, although a 
Gaussian approximation might be suitable. 

We have studied the distance d as a function of g when Q = 0.49 in the two-vortex, four-vortex, and seven-vortex 
systems studied in §3. Figure 9 displays the distances between the nearest-neighboring vortices in systems with 
N = 2,4, and 7. The vortex patterns are stable only in small parameter regions with small values of g. The stable 
parameter region for N = 2 is particularly small. Figure 9 suggests that the vortex number N increases with g. When 
g is increased beyond the critical value on the right, the vortex nucleation occurs from the surface- mode instability 
and the vortex number increases. When g is decreased beyond the critical value on the left, vortices disappear as a 
result of the symmetry-breaking instabilities. The vortex patterns with different N are multi stable in overlapping 
parameter regions. These behaviors are similar to the case of g = 1 studied in §3. However, we have not succeeded 
in evaluating the distance between the nearest-neighboring vortices theoretically, because a suitable ansatz for the 
wave function other than eqs. (15) or (16) has not been found yet. Note that the g dependence of distance d does not 
appear in the point- vortex theory 



9 



3.5 



-a 




2 J . . . . 1 

0.01 0.02 0.03 0.04 0.05 

g 

FIG. 9: Distances between the nearest-neighboring vortices in systems with N = 2,4, and 7 as a function of g for = 0.49. 

VI. SUMMARY AND DISCUSSION 

We have proposed a type of Ginzburg-Landau equation. Using the model equation, we can efficiently find stable 
vortex patterns in rotating Bose-Einstein condensates. We have found that multiple- vortex patterns are stable at 
the same frequency f2 for g = 1. We have also found that symmetry-breaking instabilities occur for vortex patterns, 
which lead the annihilation of vortices. 

The symmetry-breaking instabilities are qualitatively consistent with those obtained by the point-vortex theory, 
but the critical values are considerably different. It is left to a future study to estimate the critical values by taking 
the effects of the harmonic potential and the finite core size into consideration. The forces from the mirror images 
are essential for the symmetry-breaking instabilities in the point-vortex theory. From the analogy with the point- 
vortex system, we think that the confinement effect due to the harmonic potential is very important for causing the 
symmetry-breaking instabilities in the Ginzburg-Landau type equation. 

We have verified that symmetry-breaking instabilities also occur in a type of complex Ginzburg-Landau equation 
when the energy dissipation is small. Spiral-like trajectories have been numerically obtained, although the trajectories 
are rather complicated. We expect that such spiral-like trajectories of vortices might be observed in real experiments. 

The quantitative agreement between the point-vortex theory and the numerical simulation at g = 1 was not 
satisfactory. In the limit of large g, the core size of a vortex becomes infinitesimal, because the core size is scaled 
as 1/iJg for large g from an argument based on the scale transformation. This is intuitively because the stronger 
repulsive interaction for larger g makes the density \4>\ 2 more uniform and therefore the core size smaller. The point- 
vortex approximation is therefore expected to improve in the case of large g. We have also performed a numerical 
simulation near the L-L-L regime, where ~ Vk and g is sufficiently small in contrast to the limit of the point- vortex 
approximation when g — > +oo. We have found stable vortex patterns with N = 2,4, and 7, when g is sufficiently 
small. The single vortex is well approximated by eq. (14); however, the ansatz eq. (16) is not satisfactory for a 
two- vortex system. It is left to future studies to find a better approximation of (f> in this parameter range and obtain 
a theoretical value for the distances among vortices. 
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